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Abstract 

A tool for computing a jitter/stability metric used in NASA re- 
quirements statements is developed. An efficient algorithm is given 
for computing this metric. Two ways of implementing it on a comput- 
er are discussed. One is optimized for computational speed while the 
other sacrifices some speed to conserve memory. Timing studies are 
given to show that the improvement of computation times using the 
present algorithm over previously existing techniques can run to sev- 
eral orders of magnitude, and that previous techniques were so costly 
that the present algorithm represents enabling technology. Further 
comparisons show that the memory conservative implementation runs 
at about half the speed of the fast implementation, but can cut the 
major data storage requirement of the fast implementation by 95-99%, 
making the algorithm implementable on much smaller computers, such 
as PC’s, than it would be otherwise. Software for both implementa- 
tions is included in version 2 of the NASA time and frequency domain 
analysis program PLATSIM. 
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Nomenclature 


a, b , c, d 
J(y,w) 

A', Ax, • • • , k m 
l 

m 

m i 


n, rii 
N 

Al , ^2 , • • • , i-n 


T 

til, . . . , u m 
W, Wi,..., w m 

y 


y(i) 

y(U) 


Constants used in the jitter calculation timing formula, sec 
See Definition 1 

Number of data points in a time series sampled at equal time increments 
which are covered by windows uq, . . . , w m 

The number of different time signals being simultaneously analyzed by 
the memory conservative implementation of the algorithm 
Dominant minimum pointers, same units as t\ . etc. 

Number of windows for which jitter/stability is to be determined 
Number of windows for which jitter/stability is to be determined in 
timing test case number i 
Number of data points in a time series 

An integer parameter used by the memory conservative implementation 
of the algorithm in managing dynamic memory allocation 
Sample times for a discrete time series; arbitrary time units can be 
used, for a spacecraft instrument boresight pointing error application, 
typical unit would be sec 

Total duration of discrete time signal, i.e. , t n — t\. in the same units as 
t\ . etc. 

Dominant maximum pointers, same units as t\ . etc. 

Lengths of jitter/stablility defining windows, same units as t\, etc. 
Discrete time series of n points; units are problem dependent, for a 
spacecraft instrument boresight pointing error measurement, typical 
unit would be arcseconds 

Same as y{U), this notation only used if time points are equally spaced 
Element i of time series y , same units as y 
See Definition 1 
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( w ) 

Vl 

( w ) 

Vu 

Zi 

Greek 

n 

Subscripts 

i 

i 

j 


= See Definition 1 
= See Definition 1 

= Calculation time per data point of test case /, sec 

= Execution time of test case number i, sec 

= Time point number, 1 < i < n 
= Test case number, 1 < i < 3808 
= Window number, 1 < j < m 


Introduction 

The purpose of this paper is to provide a computational tool whose 
uses include design and analysis of spacecraft. It presents an efficient method 
for calculating a metric that has been used by NASA to quantify requirements 
for pointing jitter and stability in spacecraft instrumentation [2], This met- 
ric has been used to state requirements for the Upper Atmosphere Research 
Satellite (UARS), the Earth Observing Satellite (EOS), and individual in- 
struments on these satellites. This jitter/stability metric has the advantage 
of having a high intuitive content; it is easy for an engineer to picture the 
connection between the jitter/stability value this metric assigns to any given 
time signal and the level of noise or drift in the signal. Since NASA states 
requirements in terms of this metric, it is important to be able to calculate 
it. However, anecdotal information ([3]) indicated that when attempts were 
made to calculate the jitter or stability of lengthy time signals, the compu- 
tational burden was unacceptable. The computation was found to require 
anything from overnight or week-end computer runs in order to complete 
a single analysis to computations, such as the estimated 127 day example 
presented later, which were too computationally intensive to be feasible. 

Another feature of this jitter metric is that its very definition im- 
mediately leads one to an easy algorithm for its computation. Anecdotal 
information [3] indicates that this was the algorithm in use when unaccept- 
able computational burdens were encountered. The problem is that this 
algorithm is not the most efficient way to calculate this jitter metric. The 
exact degree of inefficiency depends on the signal length, the sampling fre- 
quency, and another problem parameter (the window length); but problems 
of engineering interest have been found for which the easy algorithm is three 
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or four orders of magnitude slower than the algorithm to be presented in this 
paper. The penalty paid for the additional computational speed of the new 
algorithm is a vast increase in algorithmic complexity. 

Examples were this metric has been used can be found in industry 
working papers ([4, 5]), in conference proceedings ([6, 7]); and an occasional 
mention can even be found in the journal literature ([8]). Reference [6] says, 
“Jitter . . . requirements [are] expressed as a worst-case change in pointing 
. . . across each time interval of interest . . . .” In reference [5] “Pointing Jitter” 
is defined to be “the peak-to-peak variation of the actual pointing direction 
over relatively short time intervals” with “Pointing Stability” using the same 
definition except that the time intervals are “relatively long.” 

In this paper, the “time intervals” used as a parameter in the defini- 
tion of the jitter/stability metric will be called windows. Then the jitter or 
stability value of a given time signal with respect to one of these windows 
is found by placing the window over the time signal in such a position as to 
maximize the peak-to-peak variation of that portion of the signal which is 
under the window. It is this maximum peak-to-peak variation which is taken 
as the value of the metric. Since the same algorithm is used to calculate 
either the jitter metric or the stability metric, the word “jitter” will be used 
for the remainder of this paper to refer to either jitter or stability as defined 
above. 

The technique presented in this paper can be applied to an arbitrary 
discrete function of finitely many time points. The units will be determined 
by the specific use which is being made of it. For quantifying the jitter in 
the boresight pointing error of an observational instrument on a spacecraft, 
typical units for time might be seconds while units for the boresight error 
might be arcseconds. On the other hand, if one is looking at irregularities in 
tectonic plate drift, the units might well be centuries and centimeters. So, 
for the most part in this paper, the units will be left unstated. 

Fig. 1 shows one such window in four of its possible positions on a 
discretely sampled function of time. The length (duration) of each time 
window, represented by the width of the rectangle, is the same for all four 
window positions. The top and bottom of each rectangle is positioned to show 
the peak-to-peak variations of the function within the window. The tallest 
rectangle represents one window position where the peak-to-peak variation 
reaches its worst case value, so its height represents the jitter in this time 
function for this window size. Note that the leftmost window extends into the 
negative part of the time axis where there are no data points from the discrete 
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time signal. Including such windows does not change the value calculated 
for the jitter metric, and will be useful in implementing the algorithm. 

In order to calculate this jitter metric with respect to a given window, 
a position must be found for this window which maximizes the peak-to-peak 
variation of the time signal under the window. The easy algorithm to ac- 
complish this calculation is to simply try all window positions and pick the 
worst case from among them. If the signal is discrete and of finite length 
(discrete signals of finite length can arise if the data come from sampling the 
output of some sensor or if the data are generated by a numerical simulation), 
then there are only a finite number of window positions with the property 
that any two of them cover different sets of points; so this is an algorithm 
which can be implemented on a computer. This search can be performed in 
a systematic fashion by first placing the window at the beginning of the time 
signal and calculating the peak-to-peak variation under it, then moving the 
window forward in time until it covers one additional point of the time signal 
and repeating the peak-to-peak variation determination, and continuing to 
advance the window one point at a time and determine the new peak-to-peak 
variation until the end of the signal is reached. The measure of jitter is then 
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Figure 1: A jitter window in four positions 
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the worst of these peak-to-peak variations. This is easily programmed in 
FORTRAN using a nested pair of DO loops or in MATLAB (MATLAB is a 
registered trademark of The Matli Works, Inc.) using a vectorized determi- 
nation of the peak-to-peak variation under the window in a given position 
inside a single for loop. 

If the data points in a discrete signal are equally spaced, the compu- 
tational complexity of the jitter calculation for each window is determined 
by two parameters: n, the number of points in the data sequence; and k, the 
number of points of data that lie under the chosen window. These numbers 
can be quite large. For example, the data may come from a simulation driven 
by disturbances with high-frequency components, so that accuracy requires 
an integration interval which is small compared to the window size, resulting- 
in a large value for k. A small integration interval will also drive up the value 
of n, as will a need to consider long term effects. 

If the calculation of jitter is programmed on a computer using the 
easy algorithm, the time to calculate the peak-to-peak variation of the signal 
under the window in one position is O(k), and n— 7 + 1 window positions must 
be considered, so the total calculation time is 0(k(n — k + 1)). If n is large 
and 1 C k <C n (in one jitter analysis of the EOS AM-1 spacecraft [7], values 
of n and k on the order of 10 5 and 10 4 , respectively, were common), this takes 
a substantial amount of time. In an analytic jitter study, the computational 
expense can be driven up by such factors as: (a) Each output signal may 
need to be analyzed over several different windows (4-7 in the EOS AM-1 
study); (b) For each disturbance scenario, many output signals need undergo 
the Step (a) analysis (20-30 in the EOS AM-1 study; representative times 
for a jitter calculation at this level using both the easy algorithm and the 
one given in this paper are given in Table 1); (c) Step (b) must be repeated 
for many disturbance scenarios (over 20 in the EOS AM-1 study); (d) Step 
(c) must be repeated for each design iteration. At this point, the need for 
a more efficient jitter calculation should be apparent. At the same time, 
approximation techniques should be avoided if possible. 

The purpose of this paper is to present a much more efficient algorithm 
which exactly calculates the NASA jitter metric. This algorithm can be used 
to simultaneously calculate jitter in a time series for each of several windows. 
A list is maintained of potentially significant maxima and another of minima 
as the time series is scanned. At each time step, information from these lists 
is used to update variables which accumulates the metric of jitter for each 
window, leading, at the end of the scan, to a determination of jitter values 
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for all windows. If there are m windows, the time to do this is (am + b)n 
where the constants a and b are independent of the actual window sizes. The 
calculation time for the first window is ( a+b)n . Comparing this to the timing- 
estimate 0(k(n — k + 1)) for the easy algorithm, it is seen that the removal of 
the dependency on the window length k has replaced a timing formula which 
is jointly quadratic in n and k by one which is linear in n and independent of 
k. The incremental calculation time for each additional window is an. This 
incremental time for each additional window has been empirically observed 
to be about 30% of the time for the first window, so the computational time 
to calculate jitter for several windows simultaneously is substantially better 
than doing the calculation for each window individually. 

This algorithm has been implemented in software and included in 
version 2 of the PLATSIM software package ([9]). PLATSIM is a MATLAB 
based software package, developed at NASA Langley Research Center, which 
performs time and frequency domain analysis of the response of a controlled 
or uncontrolled flexible structure to disturbances. Version 1.0 of PLATSIM 
([10], [11]) contains an earlier version of this algorithm (without the memory 
conservative option). 

First, the algorithm is given. Then, notes on two implementations of 
the algorithm are presented. The first implementation emphasizes computa- 
tional speed while the second implementation sacrifices some computational 
speed in order to conserve computer memory. Results of timing studies to 
quantify the improvements in computational time and to show the depen- 
dency of timing on problem parameters are shown. A few remarks on the 
amount of memory saved using the memory conservative implementation of 
the algorithm are made. The paper ends with a section of conclusions. 

The jitter calculation algorithm 

In this section the problem is dealt with of calculating the jitter in a 
time series over one or more jitter windows. Although the windows under 
which the maximum peak-to-peak excursions of the time series are being- 
observed are thought of as sliding continuously down the time axis; for pur- 
poses of determining jitter, a window need only be considered when it is 
positioned so that its right end is coincident with one of the points of the 
time series. (A window in this position is referred to as being in standard 
position ; the windows shown in Fig. 1 are placed in standard position.) This 
can be seen by observing that if a window in standard position starts sliding 
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to the right, the maximum peak-to-peak excursion under that window can- 
not increase until a new point of the time series comes into the window. This 
happens exactly when the window reaches the next standard position. So, no 
jitter information is lost by ignoring the window in any of the intermediate 
positions. 

Notation Let n be the number of points in the time series and denote 
the time series by {y{U ) |1 < i < n}. Assume that t\ < tv < ■ ■ ■ < t n - 
Suppose that there are m windows and they have lengths W\ < ••• < w m . 
If the time points are equally spaced, the time series might be ivritten as 
{y( 1), t/(2), . . . , y(n)} and kj used to denote the number of points covered by 
a window of length Wj in standard position. 

In Fig. 1, a time series is shown with n = 51, t\ = 0, and t n = 3. One 
window is shown (in four different positions), so m = 1. For this window, 
W\ = .48, and so k\ = 9. 

Definition 1 For a window of length w, denote the running tallies of the 
upper and lower bounds of y wider each window position by and y ^ 

and define them by 

VuHti) = max{y{tj)\ti - w < tj < L} (1) 

y ( L ] {ti) = mm{y{tj)\ti - w < tj < U} (2) 

where 

1 < i < n. 

Denote the jitter time history of y for this window by y ^ 

it by 

y ( j w) {ti ) = y ( u ] {ti) - yPtfi), i< i <n 
Finally, denote the jitter by J(y, w), and define it by 

J{y, w) = max{y ( j w) {ti)\l < i < n} (5) 

Where the window length is clear from context, the superscript (iv) 
will be omitted. 


(3) 

and define 


(4) 
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In Fig. 1, the windows shown correspond to i = 7, 27, 39, and 51. 
The numbers yu{U ), yhiU)-, and yj(U) represent, respectively, the top y- 
coordinate, the bottom v-coordinate, and the height of the corresponding 
rectangle. Note that even though k = 9, the leftmost pictured window covers 
only 7 points; its left side extends into negative time where the time series is 
not defined. 

One easily sees that, in Eq. (5), the same value will be computed for 
J(y. w ) if the lower limit of 1 is changed to a larger number, so long as t\ 
lies under the first window used. If jitter were to be calculated using the 
easy algorithm, this lower limit would be A;, the number of points covered 
by a window of length w positioned at the beginning of the time series. 
Doing the calculation by using the easy algorithm refers to actually doing 
the calculations in equations Eqs. (1), (2), (4), and (5). Referring to Fig. 1, if 
one wished to calculate jitter using the easy algorithm, it would be inefficient 
to use the first window position shown there. It would be more efficient to 
start with a window positioned with its left edge at time 0. If the time steps 
are evenly spaced so that there are k points under a window in standard 
position, then doing the calculations in equations Eqs. (1) and (2) requires A: 
references into the y vector for each entry in y L r and y L . Since the first A: — 1 
entries of yj are dominated by |/j(A:), these need not be computed, and only 
the last n — k + 1 values of yu and yL are needed. 

The comparative speed of the algorithm to be presented in this paper 
can now be explained. When jitter is calculated using the easy algorithm, 
then for typical values of n and A:, most entries in y are being referenced A: 
times. The algorithm presented here reduces the number of references per 
time series entry from k references (for most time series elements) to a small 
number of references which is bounded independent of A:. Some overhead is 
incurred on each pass through a time series. This is amortized over all the 
windows being considered in the pass, so the bound on the number of times 
a given series element is referenced per window improves as more windows 
are considered simultaneously. 

If one only wishes to calculate J(y,w ), it is not necessary to spend 
computer memory to contain the full arrays yu , 'IJl, and yj. For each i, 
one could calculate the 7-tli element of each of these arrays, accumulate the 
necessary information from this calculation in the jitter-to-date tally, and 
discard these i-\ li elements. 

The following discussion details the calculation of //; (/;!• The calcu- 
lation of yL is parallel. One need only reverse the direction of inequalities in- 
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volving elements of y and replace “decreasing” by “increasing” in statements 
about the monotonicity of potentially significant extrema lists of elements of 

y- 

The outer loop of the easy algorithm can be thought of as making 
a pass through the time signal. Each step of this pass places the scanning 
window in a new position. The inner loop then scans every element of the 
time signal under the window to determine the maximum and minimum 
so the peak-to-peak variation can be computed. If this calculation needs 
to be done for more than one window, the calculation for each window is 
traditionally done independently of that for the others. 

The present algorithm also makes a pass through the y vector. Howev- 
er, instead of scanning the length of a window back from the current position 
to determine the maximum and minimum under that window, the informa- 
tion needed to find that maximum and minimum has already been recorded in 
some bookkeeping arrays in a manner which allows for efficient retrieval and 
update. Another advantage to this technique is that the jitter calculation can 
be done for several windows simultaneously at less computational cost per 
window than performing the calculation for each window individually. Two 
lists are being maintained, one a list of potentially significant maxima and 
the other a list of potentially significant minima. The potentially significant 
maxima (respectively, minima) list contains elements of the y vector which 
have already been scanned and which have the potential of being the maxi- 
mum (respectively, minimum) element in one of the windows in the present 
and/or some future position of the pass. Each window has three bookkeeping- 
items associated with it, a jitter-to-date tally, a dominant maximum pointer, 
and a dominant minimum pointer. The jitter-to-date tally records how much 
jitter there is in the time series up to the point most recently scanned. The 
dominant maximum pointer points at the location in the list of potentially 
significant maxima which gives the maximum value of the time series under 
the window when its right end is at the current element of the pass. There 
is an analogous relationship between the dominant minimum pointer and 
the list of potentially significant minima. All of these bookkeeping elements 
must be updated with each step through the y vector. The efficiency of the 
algorithm comes in that this update can be accomplished in an average time 
bounded by an expression of the form am + b where a and b are absolute con- 
stants which are, in particular, independent of the lengths of the individual 
windows. 

The process of making one step through the y vector can be thought of 
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as occurring in the following sequence: Initially, all windows have their right 
endpoints at the last element of the y vector which has been scanned; for each 
window, the jitter in y up to the last scanned point has been entered into a 
jitter-to-date tally: all the information in the potentially significant extrema 
lists is relevant to windows in this position and future positions; and the 
dominant element pointers are correct for this window position. Then, the 
next point of the y vector is scanned and the windows all slide from their 
old position to their new position with their right endpoints all coincident 
with the time at which the new y point occurs. Then, the new significant 
local extrema lists are formed from the old. Then, the dominant element 
pointers are updated to be correct for the new window position. Finally, the 
jitter-to-date tally for each window is updated with new information from 
this step. 
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Figure 2: Potentially significant maxima and dominant element pointers be- 
fore update 

The process of advancing the calculation of yu by one time step can be 
observed from Figs. 2 and 3. In both figures, the same generic time sequence 
y is shown. It is sampled at 17 equally spaced points, so n = 17. There are 
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four windows with respect to which jitter is being calculated (m=4). The 
window lengths are given by the point counts Ay = 3, A: 2 = 7, A ’3 = 11 , and 
hi = 15. Fig. 2 shows the situation in computing y v after the scan of y has 
reached y( 15) while Fig. 3 shows the situation after reaching y( 16). The four 
windows overlap, having a common right edge at the dashed line labeled R. 
The left edge of w\ is at the dashed line labeled LI, the left edge of W 2 is at 
L2, and so on. Potentially significant maxima are shown as stars, other data 
points as x’s. Dominant element pointers are shown for each window. Thus, 
Vu l] ( 15) = 2/(14), ^(15) = y{ 9), etc. 

At the point of the calculation shown in Fig. 2, the algorithm has only 
seen y( 1), . . . , y( 15) and it is computationally efficient for the algorithm to 
assume that the series goes on indefinitely. One possibility the algorithm 
must consider is that the series is monotone decreasing for at least another 
14 points. If that is the case, every one of the potentially significant maxima is 
either the maximum element in one of these windows in the present position 
(these are marked by the dominant element pointers) or will become the 
maximum element under one or more windows at some future position (or 
both). The possibility that it might become a maximum in some future 
window is the potential significance of each of these elements in computing the 
maxima of the series y over some jitter windows. Note that the subsequence 
of y formed by the potentially significant maxima is monotonically decreasing 
and that each dominant element pointer points to the left-most potentially 
significant element under its window. This always happens, and can be 
established by induction on i. 

Now observe what happens when the scan of y advances to i = 16 in 
Fig. 3. The newly scanned element always becomes potentially significant. 
In this example, it is bigger than several of the potentially significant max- 
ima at the previous stage (located at i = 15,14,12, and 11), so they lose 
their potential significance. Note that the demoted elements were placed 
consecutively at the right end of the old list of potentially significant max- 
ima, This always happens; and because of this, the demotion of previously 
potentially significant extrema is an efficient computation. The dominant 
element pointers are updated. For window 1, the previous dominant element 
lost significance, so the pointer was redirected to the newly scanned element 
2/(16). For windows 2 and 4, the previous dominant element dropped out of 
the left end of the window as it moved on, so the pointer was moved to the 
right in the list of potentially significant maxima until an element was found 
in the window. Since points are equally spaced, the pointer needed only to 
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move to the next element in the list of potentially significant maxima. For 
unequally spaced points, several steps might have been necessary. The third 
window could leave its dominant element pointer where it was the previous 
step. Note that at this point, the same element, 7/(7) , is the maximum under 
both the third and fourth windows. 

The reader should note that in one more step ( i = 17), the only po- 
tentially significant maximum will be the new point 7/(17), and all dominant 
element pointers will point to it. 

The efficiency of the algorithm comes from the way in which the lists 
of potentially significant extrema and the dominant element pointers are 
updated from step to step. As another step is taken, the new element from 
the y vector is potentially significant as both a maximum and minimum. 
Whether either of these potentials is realized will depend on future behavior 
of the y vector. So the new element will be added to the end of both lists. 
What is important is that the presence of this new element on a list may 
remove any potential significance which some of the previous elements had, 
so they are deleted from the list. 

Specifically, when a new element is added to the list of potentially 
significant maxima, then any element already on the list which is less than 
or equal to the added element is no longer potentially significant. This loss of 
potential significance is seen by observing that a window in its new position 
or in any future position which covers the earlier, dominated, element also 
covers the element just added, so the earlier element is not needed to calculate 
the peak-to-peak excursion in any such window. Thus, the list of potentially 
significant maxima will be updated by removing all elements less than or 
equal to the new element and then adding the new element at the end of 
the resulting list. Consequences of updating the list in this way are that the 
elements in this list form a monotonically decreasing sequence, and that the 
set of elements deleted from the list as the result of the new addition are 
located consecutively at the end of the pre-update list. These consequences 
may be demonstrated by induction on the position of the scan in the y vector 
(at the initial step, the list contains a single element and so is vacuously 
monotonic decreasing). The list of potentially significant minima is updated 
similarly. 

Since the list of potentially significant maxima is a decreasing se- 
quence, the dominant maximum pointer for each window points as far to 
the left in the list as it can be and still have the element it points to lie in 
the window in its present position. This characteristic of the position of the 
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dominant maxima simplifies the process of updating their pointers after each 
step. 

The dominant element pointers are updated next. For each pointer, 
there are three possibilities. The first possibility is that the element being 
pointed to by a given window’s pointer may have been removed because it is 
dominated by the new element. Then the pointer is redirected to this new 
element. The second possibility is that the element being pointed to by a 
given window’s pointer may no longer be in that window after the window 
has moved to its new position. This possibility happens when, before the 
update, the dominant element was at the left end of the window. In this 
case, the potentially significant element list is scanned starting at the old 
pointer position and moving to the right until an element is found which is 
in the window in its new position. If the time points are equally spaced, it 
is only necessary to move one element down the list. Third, if neither of the 
two previous circumstances hold, the pointer is left where it was. 

For a given window, the maximum and minimum of y over the window 
in its new position are determined by the elements in the potentially signif- 
icant extrema lists pointed to by the window’s dominant element pointers. 
From these values, the peak-to-peak excursion of y over the window in this 
position can be determined. This is compared to the jitter-to-date tally for 
this window, and if the new peak-to-peak excursion is greater, the tally is 
updated with this value. 

Once the step involving the last point of y has been completed, the 
desired jitter values are found in the jitter-to-date tallies. This informal 
algorithmic description is formalized in algorithmic language as follows: 

Algorithm 1 Step 0. Initialization. Initialize the potentially significant 
maxima and minima lists to {y{t-i)}. For j = 1, 2, • • • , rn, set y^^ifii) = 
y^ 3 \t i) = y{t\), and set y'j’ 3 \t\) = J(y,Wj ) = 0. Set the iteration 
counter, i, to 1 and set all the dominant maximum pointers u\, ■ ■ ■ , u m 
and all the dominant minimum pointers l\, ■ ■ ■ , l m to t\. 

Step 1. Increment i by 1. If i > n, terminate; the desired jitter values are 
in J(y,Wi), ...,J(y,w m ). Otherwise, continue. 

Step 2. Start scanning the potentially sigiiiftcajit maxima list from right to 
left. Each time a scanned element is less than or equal to yltfi, re- 
move it. Stop the scan the first time an element is encountered which 


15 



is greater than y(tj ) or when the list becomes empty. Perform the same 
operation on the potentially significant minima list reversing the in- 
equalities. 

Step 3. Add y(tj ) to the right end of both potentially significant extrema 
lists. 

Step 4. For j = 1, 2, • • • , rn; 

Step 4.1 Update the pointers to the dominant maxima and calculate 

(a) if y (u j) was removed from the potentially significant maxima 
list at Step 2, set Uj = fi; otherwise, 

(b) if Uj < U — Wj, then scan the potentially significant maxima 
list starting at y(uj ) and moving to the right until the first 
element y(t s ) is found for which t s > U — Wj, and set Uj = t s ; 
otherwise, 

(c) Uj remains the same. 

Then set y^^ftf) = y{uj). 

Step 4.2 Similarly, update the pointers to the dominant minima and 
calculate y^ 3 \u)- 

Step 5. For j = 1, 2, ■ ■ ■ , m, set y { p\u) = yu^itfi-y^iU) and ify { p\ti) > 
J{y,Wj), then set J{y,Wj ) = y ( p ] {ti). 

Step 6. Return to step 1. 

Note: If time points are equally spaced, then the bookkeeping in this 
algorithm can be reduced by using the subscript, t. instead of the time value, 
U, for the dominant element pointers and using the number of points, Ay , 
under a window instead of its length in seconds, Wj. in the Step 4.1(b) and 
4.2(b) tests. 

Starting the algorithm at i = 1 would be wasteful if jitter were being 
calculated using the easy algorithm. However, in this algorithm, it serves 
the very useful purpose of properly initializing the potentially significant 
extrema lists and dominant element pointers, and it simplifies the additional 
bookkeeping necessary for tracking multiple windows simultaneously. 
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Remarks on cost: For the most part, these remarks address the 
case of equal time increments. The time required to calculate jitter using 
the easy algorithm for a single window covering k points would appear to 
be proportional to number of times an element of the y array is referenced. 
Thus, the time to calculate jitter using the easy algorithm is proportional 
to k(n — k + 1). If the time series y has duration T (e.g., if y( 1) occurs at 
time 0 and y(n) occurs at time T) and the window has length (duration) 
w. then k is the smallest integer such that k > (n — 1 )w/T. So, for fixed 
values of T and w, the time to calculate jitter using the easy algorithm is a 
quadratic function of the number of points, n, used in discretizing the time 
interval. The coefficients of this quadratic function depend on the length, 
w. of the window. The total time needed to calculate jitter with respect to 
several windows is found by adding up the times to calculate jitter for each 
individual window. 

On the other hand, an examination of Algorithm 1 shows that the time 
to calculate jitter using this algorithm may be bounded by an expression of 
the form amn + bn + cm + d where n, 6, c, and d are absolute constants. The 
cm + d part of this comes from Step 0 of the algorithm and, in any realistic 
example, is an inconsequential part of the total time. The total time for Step 
1 is proportional to n. The time in Step 2 is spent either in removing elements 
from lists of potentially significant extrema or in examining an element of 
one of the lists and discovering that it stops further removal. Since a given 
element of y can be removed from a list only once during the entire course 
of the jitter calculation, the “remove” part takes time proportional to n. 
Since it only takes one element per step to stop a removal back scan, this 
part of Step 2 also takes time proportional to n as does Step 3. Step 4 is 
executed n — 1 times, so Step 4.1 executed is m(n — 1) times. Most of the 
parts of Step 4.1 require an execution time which can be bounded by a fixed 
amount. The only exception is the part of substep (b) referring to a scan 
of elements in a potentially significant extrema list, and that only when the 
y points are unequally separated in time. But here, considering the entire 
jitter calculation, a given element of y can only be scanned once per window 
per list. Thus, the total time spent on Step 4.1 is O(mn). Exactly the same 
considerations apply to Step 4.2. Similarly, the time for Step 5 is 0(mn ) and 
for Step 6 is O(n). 

Comparing the two times, it is seen that the the easy algorithm timing- 
expression is a sum of m quadratics in n, so can be thought of as being 
0(mn 2 ), while the timing expression for the present algorithm is 0(mri). For 
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large n, the savings potential is great. The extent to which that potential is 
realized will be detailed in a subsequent section. 

Implementation notes 

Two implementation scenarios are considered. In the first, it is as- 
sumed that the complete time series y is available to the jitter algorithm at 
each step in the scan. Thus, the jitter analysis code can keep track of any 
past information it needs from y by an indirect addressing scheme into the 
y vector itself. This first scenario produces the faster jitter analysis code. In 
the second, it is assumed that the points of y are fed to the jitter analysis code 
one at a time and that the jitter analysis code has complete responsibility for 
any necessary recall. The emphasis in this latter case will be on managing 
computer memory to keep memory usage within reasonable bounds. 

While programming for speed of execution hardly needs justification, 
in this era of multi-megabyte computer memories the need for care in memory 
usage may. Consider, then, a time simulation involving a 1 kilohertz sample 
rate extending over a 1000 second duration. Suppose that 50 different out- 
puts must be analyzed for jitter. These numbers are within the bounds of 
parameters that have arisen in actual applications of this technology. Then, 
each output time history contains 1,000,001 double precision numbers, which 
will be assumed to occupy the IEEE standard of 8 bytes of computer storage. 
Taking into account that all 50 output time histories must be generated and 
stored at one time (running multiple simulations to generate the output time 
histories in smaller groups has its own penalty in terms of execution time) 
it is then seen that just storing these 50 y ' s requires over 381 megabytes of 
computer memory (recall that computer manufacturers use “megabyte” to 
refer to 2 20 bytes, not 10 6 ). While it might be possible to equip a high-end 
workstation to accommodate this in physical memory, a much more common 
result of trying to run a problem this size would be that a lot of virtual 
memory paging would be necessary. This would cause a severe increase in 
wall-clock time of execution. 

The Fast Implementation 

In this scheme, multiple outputs are analyzed one at a time so that 
the arrays used in potentially significant extrema processing can be reused 
for each output. It will be seen that this prevents what would otherwise 
be a doubling of storage needed (assuming that an integer requires half the 
storage of a double precision number) to hold just the output time histories. 
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The potentially significant extrema lists are realized by arrays of inte- 
gers which point to the appropriate elements of the y vector (i.e. , an indirect 
addressing scheme). These arrays must be the same length as the y vector. 
If y is strictly decreasing in time, then in Step 2 of Algorithm 1, no elements 
of y are ever removed from the potentially significant maxima list, so it must 
be big enough to hold them all. An increasing y imposes the same restriction 
on the list of potentially significant minima. For each list, there is an integer 
variable containing a pointer to the current active right end of the list, and 
an array of rn integers which contain pointers to the place in the list point- 
ing to the dominant elements for each window (in effect, a doubly indirect 
addressing scheme into the y vector). One of the efficiencies of the algorithm 
comes from the fact that once an entry is made into one of these potential- 
ly significant extrema pointer lists, it is never moved (although it may be 
overwritten). “Step 2” removal of dominated elements is accomplished by 
redefining the “right end” pointer to point farther to the left. “Step 3” is 
accomplished by incrementing the “right end” pointer by 1 and storing the 
current value of i in the new right end of the list. The “Step 4.1(a)” decision 
can be made by comparing the dominant element pointer to the new “right 
end” pointer and the “Step 4.1(b)” decision can be made by comparing the 
value pointed to by the dominant element pointer, the current value of 
and the number of points in the relevant window. Similar remarks apply to 
Step 4.2. The information needed to complete “Step 4” and do “Step 5” 
comes from following the dominant element pointers into the lists of point- 
ers to potentially significant extrema and following those pointers into the y 
vector. 

This scheme has been implemented in FORTRAN 77 code subrou- 
tines. Further code provides an interface to MATLAB. All this code is com- 
piled into a MATLAB MEX-file, which is a file of compiled binaries which can 
be called directly from MATLAB. In this form, it is part of the PLATSIM, 
Version 2.0, software package [9]. 

The Memory Conservative Implementation 

In order to keep memory usage to a relatively minimal level, dynamic 
memory management such as is available in the C programming language is 
employed. The scheme must be capable of calculating the jitter in several 
time series (e.g., multiple output channels) simultaneously. 

For each time series, the potentially significant maxima list and the 
potentially significant minima list will be maintained using a doubly linked 
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linear list of storage structures. By a “storage structure” is meant a data 
storage construct (such as the C programming language data type struct) 
which can hold several items or “fields” of data, possibly of dissimilar types. 
What makes a collection of these structures a “doubly linked linear list” is 
the logical organization of this set of structures into a linear ordering with 
one being thought of as the first (pictured as the leftmost) possibly followed 
by a second, a third, and so forth until a last (the rightmost) is reached. 
The programming mechanism for this logical ordering is the inclusion of two 
fields in each structure, one of which contains a pointer to the structure 
immediately to its left (or a special marker if the structure is at the left 
end of the list) and the other of which contains a pointer to the structure 
immediately to its right (or a special marker if the structure is at the right 
end of the list). Each storage structure must be able to store the dependent 
value y(ti ) of an element of the time series and also the time t, or, in the case 
of equally spaced points, the time index i at which the time series takes on 
this value. 

A first-in, first-out list of these structures will be maintained as a 
reserve stack to provide new structures as needed and on which to place 
unneeded structures as the information in them becomes obsolete. 

For each potentially significant extrema list, one pointer will keep 
track of its right end and, for each window, an additional pointer will keep 
track of the leftmost structure in each list which contains information from 
under this window (these correspond, respectively, to the “right end” pointer 
and the “dominant element” pointers of the Fast Implementation). As each 
new step in the jitter calculation is taken, each list will be updated using 
a new data point from one of the time series. Any structures containing 
information obsoleted by the new value will be moved to the input end of 
the reserve stack. Structures will be moved from the output end of the 
reserve stack to the right end of each potentially significant extrema list to 
hold the new value. The “right end” pointers will be updated to point to the 
structures just added, and the “dominant element” pointers will be updated 
as necessary. 

It is important to realize that when a structure is spoken of here as 
having been “moved” from one list to another, this does not mean that there 
has actually been a movement of data, from one part of physical computer 
memory to another. What is altered are the values of various pointers which 
define the position of the structure in one or another of the logical lists. 

If there are l time signals being jitter analyzed, then 21 structures 
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must be moved from the reserve stack at each jitter calculation step. Before 
that step is taken, a check will be made of the number of structures in the 
reserve stack. If it is fewer than 21 + 1. an attempt will be made to salvage 
structures containing unneeded information from the potentially significant 
extrema lists. If this does not bring the reserve stack up to at least 21 + 1 
structures, an additional N structures are added to the reserve stack using- 
dynamic memory allocation. The parameter N is 21 or larger so that the 
stack will contain at least 21 + 1 structures. The memory optimal practice 
would seem to be to add just enough structures to bring the total to 21 + 1. 
It is hoped that using the proposed scheme with a larger value of N will 
reduce overhead by reducing the number of dynamic storage allocation calls 
necessary. It is only suboptimal from memory usage by at most N — 1 
structures over the entire run of the jitter calculation. Reducing the number 
of dynamic storage allocation calls probably has its own storage usage reward 
as well as reducing execution time. The system must use additional memory 
to keep track of memory allocated at each call. This additional information 
would be needed when the system frees up the allocated memory. 

This stack management scheme has two important effects: (a) The 
reserve stack can never become empty. This simplifies the stack maintenance 
functions, (b) Structures removed from the lists and put back on the reserve 
stack as a result of Algorithm 1, Step 2, are not reused during the current 
jitter calculation step update process. This means that structures pointed 
to by “dominant element” pointers retain their “y” information and can be 
queried even if they have been removed from their list and placed on the 
stack. The program uses this “v” information to determine whether such 
an structure has, in fact, been removed from the list so that the “dominant 
element” pointer can be properly updated. 

The scheme to salvage unneeded structures from the potentially sig- 
nificant extrema lists is based on the following ideas. If a structure on a 
potentially significant extrema list contains information which comes from a 
point in its time series which has been bypassed by the dominant element 
pointer of the longest window, then it can be returned to the reserve stack. 
Such a point can be seen at i = 1 in Fig. 3. As the time series approaches 
its end, a structure which would need to be kept on a potentially significant 
extrema list if the time series were to continue indefinitely need not be kept 
if the information in it is dominated by that in an earlier structure which 
will never fall off the left end of its window before the end of the time series 
is reached. Again in Fig. 3, the potentially significant maxima at i = 7 and 
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8 fall into this category. Whenever additional structures are needed on the 
reserve stack, any of these structures which can be salvaged will be; and the 
need for allocating additional storage will be reassessed. 

This scheme has been realized in a C code package which integrates 
jitter analysis with the simulation code which generates the time series and 
some plot data compression code. Further code provides an interface to 
MATLAB. All of this code is compiled into a MATLAB MEX-file. In this 
form, it is part of the PLATSIM, Version 2.0, software package. 

Timing results 

Two timing studies are presented. In the first, the time to calculate 
jitter using the two implementations of Algorithm 1 are compared to each 
other and to the time to calculate jitter using the easy algorithm. These 
comparisons are made by calculating jitter in time signals of different lengths. 
In the second study, experimental verification is sought for the timing formula 
for Algorithm 1 given in the “Remarks on cost”. 

Comparison of algorithms and implementations 

The timing studies of the 2 implementations of the jitter calculation 
algorithm used data from the EOS AM-1 spacecraft [11] as packaged with 
the PLATSIM software [10]. In this study, there are 7 jitter windows (m = 7) 
of lengths 1, 1.8, 9, 55, 420, 480, and 1000 seconds. 

Two disturbances were used. The disturbance producing the smaller 
jitter analysis problem is called “MODIS scan mirror”. The time history of 
this disturbance extends over 1000 seconds at a sample rate of 200 hertz, 
so n = 200,001. The same sample rate is used to produce the 27 outputs 
to be jitter analyzed. This means that when the 7 jitter windows are in 
standard position, they cover 201, 361, 1801, 11,001, 84,001, 96,001, and 
200,001 points respectively. 

The disturbance producing the larger jitter analysis problem is called 
“MOPITT mirror scan” . The time history of this disturbance extends over 
1200 seconds at a sample rate of 1000 hertz, so n = 1,200,001. The same 
sample rate is used to produce the 27 outputs to be jitter analyzed. This 
means that when the 7 jitter windows are in standard position, they cover 
1001, 1801, 9001, 55,001, 420,001, 480,001, and 1,000,001 points, respectively. 

In order to compare execution times of alternative jitter calculations, 
Algorithm 1 was exercised in both the “Fast Implementation” and the “Mem- 
ory Conservative Implementation”, and jitter was also calculated using the 
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“Easy Algorithm” . These calculations were made in the context of a PLAT- 
SIM time domain analysis. Timing values were obtained using two timing- 
tools. In the case of the “Fast Implementation” and the “Easy Algorithm”, 
the jitter calculation could be isolated from the rest of the code at the MAT- 
LAB level, and MATLAB function cputime was used. For all cases, the du- 
ration of the total PLATSIM analysis was timed using the UNIX timing func- 
tion located, in standard implementations of UNIX, in file /usr/bin/time. 
These total times were compared to times found by running the same data 
through PLATSIM with the jitter calculation option turned off. 

The timing studies were run on a Sun SPARC 10/30 workstation. 
In an attempt to see how precise the timing numbers are, four of the cases 
involving the smaller jitter analysis problem were rerun 5 times each on a Sun 
SPARC 10/512 workstation. It should be noted that since exactly the same 
problem was run with exactly the same software 5 times, exactly the same 
amount of calculation was done each time; so the cpu time for the calculation 
should be the same for each run. However, these were both multi-tasking 
networked computers. While cputime and /usr/bin/time both purport to 
measure cpu time for the specific job, they do seem to be influenced by the 
total load on the computer. Timing runs were made overnights and week- 
ends to try to minimize the effect of other computer workload on the timing 
measurements, but some multitasking noise, probably from automatically 
operating system procedures, seems to have affected the timing results. 

The cases which were run five times each involved both implementa- 
tions of Algorithm 1, run both with jitter calculation enabled and with it 
disabled. Means and standard deviations of total run times for each case 
were calculated as were the mean and standard deviation of the isolated 
jitter calculation time in the case of the “Fast Implementation” with jitter 
enabled. From these calculations, the mean and standard deviation of the 
time to do jitter related calculations by the two implementations were calcu- 
lated. The results are given by stating the standard deviation as a per cent of 
the jitter calculation time. There was an 11% standard deviation in the time 
estimate for jitter analysis by the “Fast Implementation” found by taking 
the difference of the average time for a complete run including jitter calcu- 
lation and the average time for a complete run excluding jitter calculation. 
For the “Memory Conservative Implementation” , the standard deviation was 
41%. Most of this large standard deviation was due to a single data point 
which was well separated from the cluster formed by the other four points. 
By removing the obvious outlier from the data set, the standard deviation 
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was reduced to 14%. The standard deviation in the directly measured jitter 
analysis time for the “Fast Implementation” was 7%. 

With such a small sample size, these numbers must be considered 
to be rough approximations. However, it is probably safe to conclude from 
them that timing numbers presented subsequently are approximately correct 
in their most significant digit. 

Further softening the sharpness of timing data is their dependence 
on which compiler is used in compiling the executable modules and what 
optimization levels are called for. The results given in this paper are based on 
FORTRAN code which was compiled using the Sun “f77” compiler, version 
VI. 4, with optimization flag “-0” and C code which was compiled using 
the Sun “cc” compiler, version SC1.0, with optimization flag “-0”. However, 
when some of the “Memory Conservative Implementation” tests were rerun 
using binaries compiled using the GNU “gcc” compiler with optimization flag 
“-03”, execution times generally improved by 15-20%. 

Table 1 gives the times to calculate jitter using Algorithm 1 in both of 
its implementations and also doing the calculation using the easy algorithm. 
The “Easy Algorithm” test was made using FORTRAN 77 subroutines which 
were compiled using the same optimization flags and linked to MATLAB in 
the same way as the “Fast Implementation” test. 


Table 1: Jitter calculation execution times, seconds 


Case 

MODIS scan mirror 
(n = 200,001) 

MOPITT mirror scan 
(n = 1,200,001) 

Fast Implementation, 
direct measurement 

4.8 x 10 1 

3.0 x 10 2 

Fast Implementation, 
difference of full runs 

1.2 x 10 2 

7.5 x 10 2 

Memory Conservative Implementation, 
difference of full runs 

2.5 x 10 2 

1.9 x 10 3 

Easy Algorithm, 
direct measurement 

2.5 x 10 5 

1.1 x 10 7 fl 

Easy Algorithm, 
difference of full runs 

2.5 x 10 5 

1.1 x 10 7 fl 


a estimated 


The estimate used for the time to calculate jitter in the big problem 
using the easy algorithm was arrived at by running this case with only one 
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output signal instead of the 27 used in the rest of the cases and multiplying 
the resulting jitter calculation time by 27. This estimation was justified 
experimentally by running the small problem with 1, 2, and 3 outputs. The 
average time per output varied by about 1 part in 2700 over these three tests. 
Theoretical considerations support the hypothesis that time to calculate jitter 
using the easy algorithm is linear in the number of outputs, since the work 
to calculate jitter using the easy algorithm depends only on the signal length 
and the window sizes and is independent of the signal shape. Note that the 
estimated time to do the full problem exceeds 127 days. On the basis of this 
estimate, Algorithm 1 is claimed to be enabling technology for large jitter 
problems. 

Even considering the imprecision of timing numbers noted previously, 
some conclusions can be drawn. Note that n for the large case is almost 
exactly 6 times as large as for the small case. Times for calculating jitter for 
the large problem using Algorithm 1 are all roughly 6 times the comparable 
times for the small problem, supporting the claim that, for a given number 
of windows, the time to calculate jitter using Algorithm 1 is proportional to 
the number of points in the time series independent of the actual window 
lengths. The “Memory Conservative Implementation” seems to take about 
twice as long as the “Fast Implementation”. The differences between the 
“Fast Implementation” times measured directly and as a difference of full 
runs is at least partly explained by overhead in the PLATSIM program re- 
lated to outputting the jitter results which was not included in the direct 
time measurement; but for these two cases, this overhead should have been 
about the same. Perhaps the discrepancy is due to the imprecision and/or 
variability of timing measurements. 

The improvement over prior technology is clear. For the smaller prob- 
lem, the speed-up using Algorithm 1 exceeds 3 orders of magnitude, while for 
the larger problem the improvement in execution time is more like 4 orders 
of magnitude. 

The implications of this on turn-around time for a time domain analy- 
sis of the EOS AM-1 spacecraft using PL ATSIM are profound. Time domain 
analysis consisted of simulating the response of the 27 output signals to 
the disturbance time history, calculating the jitter in these 27 signals, and 
performing a data compression on these signals so that the PostScript files 
containing the plots of the output time histories could be made reasonably 
small. This data compression is done to realize a savings of disk space (from 
multi-megabytes per file to a few 10’s of kilobytes) and time to print the 
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hard copies (from about 30 minutes per to more like 2 minutes). The small 
problem which took almost 3 days to analyze using the old technology was 
done in less than an hour using the new. The fraction of compute time de- 
voted to jitter analysis dropped from almost 99% to about 4% for the fast 
method and about 7% for the memory conservative method (and note that 
the easy algorithm for computing jitter does not have a memory conservative 
counterpart). For the large problem, the run time decreased from an esti- 
mated time of more than 127 days to about 5 hours with jitter analysis time 
decreasing from over 99.8% of the total to about 4% for the fast method and 
about 9.5% for the memory conservative method. 

Parameter Dependency 

Theoretical analysis suggests that the time to compute jitter using 
Algorithm 1 should be of the form ( am + b)n where n is the number of points 
in the time series, m is the number of windows, and a and b are constants 
independent of all the problem parameters, specifically including the actual 
window lengths. (The “cm + <7” part of the timing formula given in the 
earlier “Remarks on cost” is in the noise level of the computer execution 
timing instrument, and will be ignored in this study.) 

A timing study was run to measure how the jitter calculation time 
using Algorithm 1 depends on n and m. and to see if it is independent 
of w i . . . w m (the individual window lengths). The study involved 3808 
individual cases. The time series length n ranged from 100,001 to 1,200,001 
in increments of 100,000. The time history which provided the input data 
was the MODIS roll output from the EOS AM-1 simulation using the same 
MOPITT mirror scan disturbance used for the larger jitter analysis problem 
in the previous subsection. For n < 1200001, the first n points were used. 
The basic window set was scaled by each of the factors 0.2, 0.6, 1.0, and 
1.3; any windows longer than the truncated time series were discarded; and 
the jitter calculation was done with each subset of the resulting window set. 
This gave values of rn ranging from 1 to 7 and 26 usable window lengths 
distributed in a fairly uniform exponential manner from 0.2 seconds to 1000 
seconds. 

These parameter dependency timing tests were made using the “Fast 
Implementation” of Algorithm 1. The test data had been previously cal- 
culated and saved for use in these tests. Timing studies were run directly 
from MATLAB and not in the context of a PLATSIM time domain analy- 
sis. Times were measured using MATLAB’s cputime function. Three timing- 
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runs were made on each of the 3808 test cases. 

The relevant data from the z-tli case are n, . the number of points in 
the time series, m,; n the number of windows, and r*, the time to calculate 
jitter. The dependent variable was taken to be z*, defined to be r, : //q, and 
the dependence of z, on m, was investigated. A scatter plot of the points 
(m,;, z, ) showed that the data tended to concentrate along a straight line, but 
a few outliers were observed. Further analysis showed that any case which 
gave an outlier in one of the three timing runs did not produce an outlier in 
either of the other two runs. It was concluded that the outliers resulted from 
external influences on the timing instrument. This was taken as justification 
that the outliers could be eliminated. 

Outliers were eliminated by taking z, to be the median value of the 
ratios formed with data from each of the 3 runs. A scatter plot of 

(//>,. now concentrated along a straight line with no wild points. 



Figure 4: Linear regression fit to window size timing data 

Using linear least squares regression on this data, it was found that 
z, fa ami + b with a = 7.56 x 10“' sec and b = 1.74 x 10 -6 sec. The 
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root mean square error in this approximation was 2.73 x 10 -8 sec which was 
about 0.6% of the mean z. t value of 4.27 x 10 -6 sec. This is a tight fit which 
experimentally validates the theoretical timing formula. A scatter plot of 
the median computational time per time series point as a function of the 
number of windows is shown in Fig. 3 together with the least squares linear 
regression fit line. 

Estimating jitter computation time as (am + 6)n, it is seen that the 
time to calculate jitter with respect to one window is about (a + b)n while 
each additional window adds about an to the time. This means that the time 
to calculate jitter with respect to each window after the first is about 30% of 
the time to calculate jitter with respect to the first window. The conclusion 
to be drawn from this is that it is substantially more efficient to calculate 
jitter over several windows simultaneously than to do the calculations for 
each window individually. 

Memory Conservation Results 

In the small problem, n = 200, 001 and l = 27. Thus, it takes 200001 x 
27 x 8 bytes of memory to store the time histories of the outputs in 8 byte 
double precision words. The jitter computation also requires 2 integer arrays 
of length n to hold the information for the lists of potentially significant local 
extrema. Assuming that an integer occupies 4 bytes of storage, and recalling 
that a megabyte of memory refers to 2 20 bytes, it is seen that the major data 
storage part of the memory requirement for the “Fast Implementation” of 
doing the jitter calculation is almost 43 megabytes. 

The storage structure used in the “Memory Conservative Implemen- 
tation” must hold one double precision number (?/(£*)), one integer (?'), and 
two pointers. In the implementation of the C programming language used 
in this study (the standard C compiler cc furnished with the Sun OS 4.x), 
these individual fields occupy 8, 4, 4, and 4 bytes respectively. However, 
the system uses 24 bytes for the complete structure. One can only speculate 
about those 4 extra bytes. Perhaps it has something to do with aligning 
structures on word boundaries. 

The choice was made to allocate these structures in blocks of 10,000. 
The small problem required 4 of these blocks. Thus, the major data storage 
part of the memory requirement for the “Memory Conservative Implemen- 
tation” of doing the jitter calculation was 1000 x 24 x 4 bytes which is less 
than 1 megabyte. The memory savings here were almost 97.9%. 
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The large problem has n = 1, 200, 001 and l = 27. Repeating the pre- 
vious calculation with these numbers, it is seen that the major data storage 
part of the memory requirement for the “Fast Implementation” of doing the 
jitter calculation is over 256 megabytes. When it was run with the “Memory 
Conservative Implementation” , 9 blocks of structures were allocated, so the 
major data storage part of the memory requirement for jitter analysis in this 
run was barely over 2 megabytes. In this case, the memory savings were over 
99.2%. 

It should be noted that the memory usage for the “Fast Implemen- 
tation” is completely determined by the length and number of time signals 
being analyzed, but that the memory usage for the “Memory Conservative 
Implementation” also depends on the signal shape and the lengths of the jit- 
ter windows. An extremely noisy signal would tend to require less memory 
while a signal with long monotone stretchs would require more. So, while 
the numbers given here for the “Memory Conservative Implementation” are 
indicative of the savings, results will vary somewhat from case to case. 

Conclusions 

A fast algorithm for calculating jitter in a time signal has been pre- 
sented. When applied to large problems typical of those found in spacecraft 
design analysis, it computes the NASA metric of jitter faster than prior tech- 
nology by 3 to 4 orders of magnitude. In the case of very large problems, this 
is enabling technology. In any case, it tames the jitter analysis calculation 
so that instead of being the overwhelmingly dominant element in the total 
calculation time of a typical analysis, it is a relatively small part. 

Two implementations of the algorithm have been described. The 
faster one takes advantage of knowledge of the entire time series when this 
information is available. The memory conservative version allows the calcu- 
lation to proceed with much less memory usage, preventing excessive page 
swapping and even disk swap space overflow for particularly large problems. 
If jitter must be calculated with respect to several windows, it is much more 
efficient to calculate the jitter for all the windows simultaneously than it is 
to do the calculation for each window independently of the others. 
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